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We discuss the possibility that dark matter axions form a Bose-Einstein condensate (BEC) due to 
the gravitational self-interactions. The formation of BEC occurs in the condensed regime, where the 
transition rate between different momentum states is large compared to the energy exchanged in the 
transition. The time evolution of the quantum state occupation number of axions in the condensed 
regime is derived based on the in-in formalism. We recover the expression for the thermalization rate 
due to self-interaction of the axion field, which was obtained in the other literature. It is also found 
that the leading order contributions for interactions between axions and other species vanish, which 
implies that the axion BEC does not give any significant modifications on standard cosmological 
parameters. 
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I. INTRODUCTION 



Identifying the origin of the dark matter of the universe is one of the priorities of modern high-energy physics and 
astronomy. So far, many particle physics models of the dark matter have been proposed (see e.g. [l[ for reviews), and 
well motivated candidates are so called weakly interacting massive particles (WIMPs) and axions. Both of them possess 
suitable properties for the dark matter in that they are non-baryonic, cold, and collisionless. However, the nature 
of cosmological behavior is completely different between WIMPs and axions. WIMPs are produced from primordial 
soup of radiations, and their population is fixed when they decouple from the thermal plasma. We can interpret them 
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as a collection of classical particles, whose velocity dispersion is determined by the decoupling temperature. On the 
other hand, axions are produced non-thermally, having a small velocity dispersion compared with the temperature 
of the thermal plasma at the epoch when they are produced. Furthermore, they have huge occupation number in 
the phase space since they are bosons. Because of these properties, we interpret them as a classical field rather than 
individual particles. This classical field of axions coherently oscillating in the field space behaves as a cold matter 



These curious properties of dark matter axions motivate the possibility that axions form a Bose-Einstein condensate 
(BEC) 3]. Indeed, the energy dispersion of axions at the production time 8u ~ O(10 _13 )eV is much smaller than the 
critical temperature for BEC, T c = (7r 2 n/C(3)) 1/3 - O(100)GeV, where n is the number density of axions . However, 
this argument relies on the following assumptions. First, the particles must be bosons. Second, their number must 
be conserved. Third, they should have huge phase space density. Finally, they must be in thermal equilibrium. The 
first three conditions are obviously satisfied for dark matter axions, but the final one, whether axions thermalize or 
not, is a non-trivial issue. 

The cosmic thermalization of axions is extensively discussed in Ref. [![. Here, the thermalization means that the 
system relaxes into a state with the highest entropy by exchanging energies and momenta between particle states. 
In order to investigate such a process, some statistical mechanical treatments are required. The usual analysis using 
the Boltzmann equation cannot be applied to this problem, since highly degenerate axions are essentially fields in 
the classical limit. Such a system of axion fields does not match the assumption of the Boltzmann analysis, where 
the system is considered to be a collection of point particles. This situation defines a peculiar regime of the many 
body system, called the condensed regime, which should be distinguished from the particle kinetic regime where the 
Boltzmann analysis can be applied. In the condensed regime, the energy transfer of the scattering process is small 
compared with the scattering rate, because the interaction occurs between highly degenerate states. Thermalization 
in the condensed regime is not well understood in comparison with that of the particle kinetic regime. 

In Ref. [H, the thermalization process in the condensed regime is analyzed by describing the axion field as a set 
of quantum-mechanical oscillators (i.e. quantum operators) and deriving the evolution equation of each oscillator. 
However, in the formalism of Q the thermalization rate is estimated by comparing the "order of magnitudes" of 
quantum operators, and the occurrence of the thermalization is confirmed only by computing the quantum- mechanical 
averages of the occupation number of each oscillator numerically, which is realized for toy models with small number of 
oscillators and particles. It is impossible to realize the actual system with huge number of axions using the numerical 
scheme described in [J]. 

In this paper, we revisit the issue of the axion thermalization. The purpose of this work is to develop a robust tool 
to describe the relaxation process of highly degenerate bosonic fields. Instead of using the approach of Ref. we 
compute the expectation value (i.e. the quantum-mechanical average) of the occupation number of axions, and solve 
its time evolution, which informs us of the change of the distribution function of axions. The computation is executed 
by using the technology which was originally introduced by Schwinger et al. Q, called the "in- in" formalism. This 
formalism enables us to calculate the time evolution of the expectation value of quantum operator in the systematic 
way. Furthermore, it can treat a state to which particle-like interpretation is not applied, if we use an appropriate 
representation for a state at the initial time. In our case, a coherent state is used to describe a state in the condensed 
regime. Using this formalism, we estimate the thermalization rate as the inverse of the time scale in which the 
expectation value of the occupation number changes its value. Whole things can be described in the analytic way, 
and it is not necessary to use numerical simulations. 

The discussion on thermalization of axions and formation of a BEC is not only the theoretical issue, but has a 
relevance to observations. There are several observational evidences indicating that the phase space structure of 
galactic halos is consistent with the "caustic ring model" @ . In this model, the high density surfaces (caustics) in the 
phase-space distribution of the dark matter particles become ring-like configurations when dark matter particles fall 
into a galactic potential with net overall rotation. Recently, it was pointed out that this caustic ring model is predicted 
if dark matter axions form a BEC 0], which might be regarded as an evidence of the axion dark matter. However, 
it was also suggested that axion BEC might enter into thermal contact with photons, and modify some cosmological 
parameters from the standard values (J,l8|- In particular, if axions and photons reach thermal equilibrium, the photon 
temperature is cooled, which predicts a smaller value of the baryon to photon ratio at the big bang nucleosynthesis 
and larger value of the effective number of neutrino species N e ff. The predicted value is N e s — 6.77, which is larger 
than the observed value N c g ~ 3-5 Q. This result seems to be disapproval of axion BEC dark matter, but we find 



1 Here, the expression for the critical temperature T c = (ir 2 n/£ (3)) 1 / 3 used in Q is applied for relativistic particles, which might not 
be appropriate for cold axions because they are non-relativistic. Even though, the condition Sui >C T c is also satisfied if we use the 
expression for non- relativistic particles, T c = (27r/m)(n/C(3/2)) 2 / 3 ~ O(10 22 )GeV where m is the mass of the axion. 
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that such photon cooling effects are fictitious. As will be shown later, axions do not enter into thermal contact with 
photons, though they form a BEC. 

The organization of this paper is as follows. In Sec. HH we introduce the method to evaluate the time evolution 
of the occupation number of axions. The expectation value of the quantum operator for the occupation number 
is computed by using the perturbative expansion in terms of the interaction Hamiltonian of the axion field. We 
perform the calculation of the leading order terms in perturbation theory, while the second order terms are evaluated 
in Appendix [5] In Sec. IIII1 implications of the results of our analysis on cosmology are discussed. We estimate the 
relaxation rate of self-interacting degenerate axions and recover the formula for the thermalization rate of axions, which 
was obtained in Ref. We will see that the thermalization rate exceeds the expansion rate when the temperature 
of the universe becomes sufficiently low, at which axion BEC is formed. Finally, summary and conclusions are given 
in Sec. El 



II. AXION FIELD DYNAMICS 



In this section, we develop the formalism to compute the time evolution of quantum occupation number of the 
axion field. Our interest is to calculate the expectation value of a quantum operator Af p (t) at a given time t, which 
represents the number of axions occupying a quantum state labeled by three-momentum p. Such a problem can 
be dealt by using the "in- in" formalism (or Schwinger-Keldysh formalism) [a|. In cosmology, this formalism was 
applied to calculate quantum contributions to cosmological correlations [l(| [H[ ■ Following such a formalism closely, 
in this work, we calculate the time evolution of the occupation number in a systematic way by use of the perturbative 
expansion. 

The outline of this section is as follows. In Sec. Ill A[ we review the formalism described in [l(| [H| and give the 
formula to calculate the expectation value of the quantum operator. In Sec. Ill Bl the mode expansion of the field 
operator is taken and the quantum occupation number is defined. In Sec. Ill CI we discuss how to represent the 
coherently oscillating axion fields as quantum states. Using the ingredients obtained in Sees. Ill Al IIIB| and III C( we 
compute the time evolution of the occupation number due to the self-interaction of the axion field in Sec. HID] Finally, 
interactions with other particles such as baryons and photons are discussed in Sec. Ill El 



A. The in-in formalism 



Let us consider the theory with a real scalar field (the axion field) </>(x, t) in the Minkowski background. The 
Lagrangian density is given by 



C 



(1) 



H[<f>(t),w(t)] = / d 3 x(irj> - C) = H [<t>(t)Mt)] + £fr|#(t),7r(t)], 
where 7r(x, t) — <f>(x,t) is the canonical conjugate. Here, we decompose H into its free and interaction terms, 



where m is the mass of the axion and Ci is the interaction term which we specify later. The Hamiltonian of the 
system is given by 

(2) 

(3) 
(4) 

(5) 
(6) 



H [4>{t),K{t)\ = Jd 
H I [cj>{t),K{t)] = - J d?xd. 



The quantum operators satisfy the canonical commutation relations, 

[0(x,t),7r( y ,t)] ^ 3 >(x-y), 
[#x,t),0(y,t)] = [7r(x,t),7r(y,t)] =0. 

Their time evolution is given by the Heisenberg equations, 

=*[ZW),7r(t)],#x,t)] > 
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These equations can be solved formally, 

4>{t) = U- x (t,t )^(to)U(t,t ), 
7r(t) ^ C/- 1 (^t )7r(t )t/(t,io), 
for some fixed time to, where U(t,to) is given by 

d 



dt 



U{t,t ) =-iH[<P(t),Tr(t)]U(t,t Q ) and U(t ,t ) = l. 



Now we move on to the interaction picture. Let us define interaction picture fields cf) 1 and it 1 such that 

^(x,t) = *[# [^(*)V(*)],^(x,*)] , 
7T 7 (x,t) =i[H [</> I (t),n I (t)U I (x,t)] , 

with (j) 1 (to) = <fr(to) and it 1 (to) = Tt(to). Solutions of these equations are given by 

<f>\t) = ^o _1 (*»*o)^(*o)fb(t,*o), 
7T J (i) = U^^toMtoWo^to), 

where Uo(t,to) satisfies the following equation, 



By noting that 



Eqs. (0) and ((TTJ» lead to 



— Uo(t,to) = -iHo[<t> I (t),Tr I (t)]U(t,t ) and U (to,t ) = 1. 



H [<t> I (t),Tr I (t)} = Ho[<t>(t ),w(to)}, 
H[<f>(t),ir(t)] = H[<p(t ),7r(to)}, 



dt 



F(t,t ) = -iHi(t)F(t,t ) and F(t ,t ) = l, 



where 

F(t,t ) = U - l (t,to)U(t,to)> 
Hj(t) = U^(tM)Himo)Mio)]Uo(tM) =H I [cf> I (t),Tr I (t)}. 
The solution of Eq. (fT3|) is given by 

rt 



F(t,t ) = Texp ^-i jf Hi(t)dtj 



and also 



^(Mo) = Texp i / H!(t)dt , 



where T (T) represents (anti-)time ordering. 

A quantum operator 0[cj)(i), ir(t)] constructed from (f> and tt can be written as 



0(t) = F- 1 (t,t )O*(t)F(t,t ) 



Texp i / 



o 7 W 



Texp -i / ffj(i)di 



(7) 
(8) 

(9) 

(10) 
(11) 

(12) 

(13) 

(14) 
(15) 

(16) 
(17) 



(18) 



where O 1 (t) = Offi (t),^ 1 (t)]. Using Eq. (fT5|) . we can compute the expectation value of the operator (0(t)) 
(\f r |0(i)|\I r ) at a given time t for an "in" state |\&) specified at the time to. It is convenient to note that 



(om = 



■N 



N=0 



to J to 



dtN I dt 



N-l ■ 



dhdHjih), [Hj(t 2 ), . . . 0'(t)] ...]]), (19) 



which can be derived by mathematical induction. 
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B. Mode expansion and occupation number 

Since the calculation in Eq. (|19|) is performed in terms of the interaction picture fields (j) 1 and it 1 , it is convenient 
to write down relevant quantities in the interaction picture. From Eqs. ([5]) and (|10p . the interaction picture fields also 
satisfy the commutation relations, 

[^(x,t) )7 r'(y,i)] = ^ 3 )(x-y), 

[^(x, t), ^(y, t)} = [tt 7 (x, t), 7r J (y, t)] = 0. (20) 
Eqs. ((9]) imply that and tt 1 are solutions of free field equations of motion, giving their mode expansions, 

=!w?^m le " "< +e '*'^' (21) 

*' <x > t) = -<I (IFFVT t e *"""p - e_ "*°?] ' (22) 



where i? p = -\/m 2 + |p| 2 , x° = t, and p° — E p . Then, the commutation relations (f2T)|) are equivalent to 

= (27r) 3 <J< 3 >(p-p'), and [a p) a p ,} = [aj? , a'J] = 0. (23) 
The creation and annihilation operators diagonalize the free Hamiltonian of interaction picture fields, 

H [<t>I(tW(t)] = J <? X [i(Tr') 2 + i(V^) 2 + imV) 2 ] = / -0^E p (a$a p + ^tt) 3 ^ (0)) . (24) 

Here, let us define the operator whose eigenvalue gives occupation number of a momentum state p, 



3 (25) 



On the other hand, its eigenstate can be obtained by applying the ladder operator a 1 ^ on the vacuum state defined 
by 

«p|0)/ = 0. (26) 

An operator similar to (|25[) in the Heisenberg picture can also be constructed. Since Heisenberg picture and 
interaction picture operators are related, 

0(x,t) = F-\t,t )<p I (x,t)F(t,t ) and tt(x, t) = F'^t, h)* 1 (x, t)F(t, t ), (27) 

the following time-dependent operators are useful, 

Op(t) =F- 1 (t,t )a I p F(t,to) and o],(t) = F- 1 (t,i )a^ J F(t,i ). (28) 

From Eqs. (|23|) and (l28l) . it is manifest that a p (t) and a£(t) also satisfy the canonical commutation relations, and 
diagonalize the free Hamiltonian of Heisenberg picture fields Ho[<f>(t), n(t)]. Hence we recognize that the operator 



A/p(i) = j^al (t)a p (t) (29) 



(2tt) 3 

describes the time evolution of the occupation number, and we substitute it into 0(t) in the left hand side of Eq. (fll 
However, in the actual calculation, the operator ([231) is used as O 1 in the right hand side of Eq. (IT^l) . 

Note that neither a p nor a p (t) diagonalizes the full Hamiltonian H = Hq + Hj, and the state |0)j given by Eq. (|2"6"]l 
is not the ground state of the full Hamiltonian, which we shall denote \0)h- On the other hand, the "in" state, which 
is used to calculate the expectation value in Eq. (fT9"|). is defined as an eigenstate of H, not Hq. Such a state can be 
constructed by applying the operator a™^, which creates one particle state from \0)h- It is assumed that this "in" 
state approaches a free particle state constructed by a^(i) in the limit to — t —> — oo, up to a factor representing the 
renormalization of the wave function. Such a factor can be absorbed into the physical mass of the field 0, which 
differs from the bare mass m appealing in Hq [l2j . 
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Aside from the renormalization factor, in the limit to — t^ — oo, the right hand side of Eq. (fl9t can be reduced into 
the expectation value in the "vacuum" |0) / multiplied by a factor arising from the overlap between states |0)j and \0)h- 
This overlapping factor drops out when we divide the expectation value by 1 =h (0|0)jj. This procedure is justified 
by taking the limit to — t — > — oo(l — ie) in a slightly imaginary direction, where e is a positive infinitesimal. Note that, 
in this case, the quantity appearing in the denominator is equal to unity because of i(0\F~ 1 (t,to)F(t,to)\0) i = 1. 
This fact implies that all vacuum fluctuations automatically vanish in the in- in formalism [lol ] . 

Having removed ambiguities via the procedure described above, we can simply calculate the right hand side of 
Eq. (1191) with the state constructed by applying operators a 1 ^ on the "vacuum" state |0)j in the interaction picture. 
This initial state will be specified in the next subsection. 

In the following, the subscript "J" is omitted for simplicity. It is convenient to consider a finite spatial box with 
volume V — L 3 so that the label of each mode becomes discrete, p n = (2ir/L)n, and n = (n x ,n y ,n z ), where n x , n y 
and n z are integers. Then we just take the following replacements, 

(27r)V 3 )(p-p') -> V5 nx , n ,J nyK 5 n „ K , 
d 3 p 1 



(2tt) 3 V 



E 



n x ,n„,n z 



and also 

[ai,Oj] = VS it j, and [a*, a,-] = [aj,aj] = 0. (30) 

Here, the italic indices i, j should be understood as abbreviated notation for three dimensional vectors with discrete 
components. The number operator defined in Eq. (|25[) is rewritten as 

Mn EE (31) 




where the factor 1/V appears due to the factor V in Eq. (|30[) . 

C. Coherent oscillation as a quantum state 

Now let us specify the state at the initial time. Since we are interested in the evolution of coherently oscillating 
axions, the initial time is set to be the epoch of the QCD phase transition, at which the mass m becomes greater than 
the Hubble parameter H and the classical axion field begins to oscillate around the minimum of the potential. These 
axions are produced due to the misalignment mechanism 0], and are called the zero mode, since a huge number of 
axions homogeneously oscillate over a large distance. 

In addition to this zero mode, however, there are additional contributions to the axion abundance. One contribution 
is produced by the thermal bath in the early universe, and its abundance is fixed at the decoupling temperature 
Another contribution comes from the decay of topological defects, such as global strings and domain walls [IT 
Both of them have definite momenta, and we call them the non-zero modes in contrast to the zero mode. If inflation 
occurred before the Peccei-Quinn (PQ) phase transition, those produced by topological defects can be a dominant 
component of dark matter. On the other hand, if inflation occurred after the PQ phase transition, their population 
is negligible. See [TH-flU for recent developments about this issue. 

Each of the zero mode and the non-zero modes corresponds to a definite quantum state. In particular, it is possible 
to construct a state with a definite momentum pk occupied by A/fc axions as a number state 

m =7mm { ° ir * loh - (32) 

where |0)j is the vacuum defined by Eq. (|26j). This is an eigenstate of the number operator pip , and we can construct 
complete orthonormal basis by using a series of the number states. Here, it should be noted that the non-zero modes 
correspond to the number states. In the classical limit, these states can be interpreted as classical point particles with 
definite energies and momenta. 

On the other hand, the zero mode has different properties compared with non-zero modes. It has a huge occupation 
number as large as Af ~ 10 61 . In the classical limit, this state should be interpreted as a classical field, rather 
than point particles [4]. Such a highly degenerate Bose gas of axions might be described as a coherent state 
Mathematically, a coherent state can be represented by using the basis of number states [2(| 



M 2 



= (4) n |0>/, (33) 
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where a, is a complex number and the numerical factor is chosen so that it is normalized (ai\ai) = 1. The coherent 
state is characterized as an eigenvector of the annihilation operator such that 

a i \a j )=V 1 / 2 a j 8 ij \a j ). (34) 

Hereafter we assume that a huge number of particles occupy a small number K of states around the ground state 
and that they are described as coherent states. There also exist non-zero modes, which occupy states with higher 
momenta and are described as number states. The collection of such states can be expressed as 

|{A0, W> = LI Tjn^ teOTW). (35) 
k>K y^k-V k 

co n 

|{a})=n^ lK|2 E^kK t )"l°^ 06) 

i<K n=0 n - V V 

Note that i < K is the abbreviated notation representing the sum over lowest K modes (i.e. actual states are labeled 
by three momenta, and we must distinguish them by spatial directions of momenta as well as their absolute value). 
Let us call the modes with i > K the particle-like modes and the modes with i < K the condensed modes. It would 
be convenient to note the following relations: 

a u (a])^] =M j V5 ij {a^-\ (37) 



(a^'Xj =Af j VS ij (a j )^- 1 , (38) 

a k \W},{*}) - | akV y l{MU a}} if k < K ' (39) 

where \{Af} k ,{a}} is the state obtained by replacing the factor (aj^* /\/NkW Nk with {a\) n "^ 1 /y/(Af k - l)!^"" 1 
in Eq. po]> . 

It is important to assume that there are plural condensed modes {K > 1). Our interest is to know how these 
condensed modes reach thermal equilibrium by exchanging their momenta. In this sense, the "zero mode" is not 
exactly a single mode with zero momentum, but the collection of K plural modes near the ground state. We summarize 
the contents of the initial state in table HI 

TABLE I: Classification of axions with their origins and quantum state representations. 

production mechanism quantum state 

zero mode misalignment mechanism coherent states (condensed modes) 

non-zero mode (topological defects) decay of defects number states (particle-like modes) 

non-zero mode (thermal axions) thermal decoupling number states (particle-like modes) 



Let us take the expectation value of <f> given by Eq. ([2T]) at the initial time to for the state {a}), 

^ = ({Af},{amx,t )\{Af},{a}) 

= -7^f(e- iEnt0+ip ^a n + e iE ^-^a* n ). (40) 

Since the wavelength of the condensed modes is comparable or greater than the QCD horizon, |p n | < H(to) ~ tjT , 
p n • x <C 1 and hence the factor e ±2Pn ' x is negligible. This approximation remains valid as long as we consider the 
dynamics inside the horizon. We also approximate E n — y/ m 2 + ~ m since the coherent oscillation begins when 
|Pn| ^ H(t ) < rn is satisfied. Then, the expectation value, 4>0i is given by 

^ £ -^( e -^o an + e imt 0<) 
n<K ^ 2mV 



51 \/-^7l"™|cos(TOt -/3 n ), (41) 



, mV 

n<K 
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with 

a n = \a n \e 1 ^. (42) 

If the condensed modes are decoupled with each other, the expectation value of the field oscillates like (<fi) oc cos(mt — 
fi n ). Each mode oscillates independently with different amplitude \a n \ and the total amplitude is given by the 
superposition of K oscillating modes. 

Next, let us take the mean square deviation of the field amplitude for a single coherent state given in Eq. ([33)) . 



-p^^k^pwf^;- (43) 

Since this result does not depend on oti, it holds for the state with on = (the vacuum state), which implies that 
the deviation given in Eq. (|43p is nothing but the vacuum fluctuation. Therefore, the coherent state has the same 
trajectory with the classical field and the same fluctuation with the vacuum. 
The expectation value of the momentum conjugate (|2"2"|) leads to 



0o= <{A0,{a}Kx,i )|{A0,{a}> 
= y'^^lanl sin(/3„ - mt ). 



(44) 



n<K 



Let us assume that the initial velocity (4>) of every mode vanishes, j3 n = mt$. In this case, we obtain 



E \/^7l a «l = E («) 



, mV 

n<K n<K 



where F a is the axion decay constant and 



2 Kl 

mV F a 



(46) 



is the initial misalignment angle for a mode n. 

The total number of particles at the initial time is given by 



iV = ^(M,WK|{Af},{a}), (47) 

TO 

where J\ n is the number operator given in Eq. pip . Dividing it by a volume V yields the number density of axions 

"tot = y = -^^({A0,{a}|a„a„|{A0,{a}) =n p + n c , (48) 

n 

where 

n>K 

is the number density of particle-like modes, and 

n c = n c,n = ynFtie™) 2 , n c , n = ^\a n \ 2 (50) 

n<K 

are the number densities of condensed modes. Here (6* ml ) 2 is the square of the total misalignment angle 

(^) 2 ^ E(o 2 - (si) 

n<K 
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In the continuous limit V — > oo, Eq. (|48|) can be rewritten as 

ntot = / (^ /(P) ' (52) 
where /(p) is the total phase space distribution function of axions, 

/(p) = A/p + J2 (2^) 3 <5 (3) (P ~ Pn)n c ,„. (53) 

D. Time evolution of quantum occupation number 

In the following, we compute the time evolution of the expectation value of the quantum number operator (fl"9]l . By 
using the in-in formalism, the time evolution of the occupation number is given by 

(#„(<)) = {tip) + i I dhdHjih),^}) - f dt 2 [ 2 dhdHjih), [fl/(t 2 ), ■#„]]) + (higher order in Hi), (54) 

where (...) represents the expectation value for the state given by Eq. (|35|) . We consider the following form of the 
interaction 

w) = ^Ei A « e " <n * , *4«, t «i«i. ( 55 ) 

where njj = E L + Ej - E k - E x and K% satisfies A% = Ag = A^ = A§*. This can be obtained from A0 4 /4! type 
interaction in the effective Lagrangian of the axion field with A ~ 0.35m 2 /F 2 , and the coefficient K k \ becomes 

Here we dropped the processes which violate axion number such as aWa^a, since the rate of such a process is too 
small and irrelevant for the case of interest Q. In addition to this self-coupling, axions also interact due to their 
gravitational potential. In the Newtonian limit, the interaction Hamiltonian of the gravitational coupling is given by 

H^W),*®] = "f / dW ^*0^ , (57) 

where G is the Newton's constant and p(x, t) — (7r 2 (x, t) + m 2 <f> 2 (x, t)) / 2 is the energy density of axions. This leads 
to the term (j55"|) with the coefficient 

A »"« - - 4 * Gm2 (r^sf + k4f) (58) 

where we used the approximation Ei ~ m. 

Let us evaluate the term of first order in Hj. Using the following relation, 

[a\a\aiaj } a) p a p } = V 5 ip a\a\aja p + Vdj p a\,a\a,ia p — V5 kp a p a]aiaj — V5i p a p a\aiaj, (59) 
the commutation relation becomes 

\Hi(t)M P ] = -L E [AXf e-^^Uk-Op - H - c ] ■ (6°) 
The following relations coming from Eq. (|3"9")l are useful to take the expectation value of the above term, 



a k ,a k \{Af},{a}) = < 



y/Mk{Mk - l)V\{Af} 2k ,{a}) if k = k>>K 

y/N k Mk'V\{Af} k - k ' ,{a}) if kj^k', k> K, and k' > K 

y/Af k ~a k ,V\{Af} k , {a}} if k > K and k' < K 

t a fc a fc ^|{7V}, {a}) if fc < if, and fc' < # 



(61) 
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where the state \{Af} 2k , {a}) contains a factor {a\) Ark 2 / 'yj '(Af k - 2)\V Mk ~ 2 for mode k, and the state \{J\f} k ' k ' , {a}) 

contains a factor {al) Ark ~ 1 {al,) Ar "'~ 1 /y/iAf k - TjW* 1 * HA4' - l)!^'" 1 for modes k and k'. By separating the 
summation over indices jkl into the contribution of particle-like modes > K and that of condensed modes < K, and 
using Eq. (IBTj) . we can compute the expectation value of Eq. (ISO")) . For p < K . after some algebra, we obtain 



c.c. 



j>K k<K 



+ 



2^ E E E [ 

j<K k<K 1<K 



AS 



kl ' 



fci a£a-; a-jO-p - C.C. 



for p < K. 



(62) 



The first term in the right hand side of Eq. (|62|) vanishes since AS contains the conservation law of three momenta 
6p+j t k+j- On the other hand, for p > K, this term exactly vanishes 



([Hi(t),fi p ]) =0 for p > K. 
Finally, taking the time integration yields the contribution at the first order perturbation, 



(63) 



2V 2 



A ki 



e -«M*( e ^S'(t-*o) _ i) 



j<K k<K 1<K 



EEE 



2^ 2 



j"<AT fc<K /<A" 



AS" r- 



for p < K, 



and 



if dt 1 ([H x (t),Ar p ]) = for p>if, 

■/to 



(64) 



(65) 



where we have dropped the rapidly oscillating term e 1 ^"^-^) as t — to —> oo in the second line of Eq. ()64[) . Note 
that, if there is no scattering (p = k or p = I), the first line of Eq. (|64|) also vanishes. 

As was conjectured in [4|, there are two distinct regimes for the interaction process. One is the particle kinetic 
regime characterized by the condition r <C 5oj, where T is the evolution rate of the system and 5ui is the typical 
energy exchanged in the interaction. In this regime we expect that Q p k \t 3> 1 and the factor exp(— iV^,t) in Eq. (|64p 
cancels out when the time average is taken. Hence the first order term (|64j) becomes irrelevant, which requires us to 
evaluate second order terms in the expansion (|54[) to follow the time evolution of the occupation number. The explicit 
calculation for second order terms is given in Appendix [X] 

The opposite regime characterized by the condition r ^> Sui is called the condensed regime. Since Qj^f <C 1 is 
satisfied, we can safely set 



(66) 



in Eq. (|64l) . and hence the first order term becomes relevant for the estimation of the evolution rate. In this regime, 
considerably small momenta 5uj are exchanged between K highly occupied states. It makes sense even though 
fffit <C 1, since the transition occurs as Afft^t ^> 1 for huge number of particles, Af. The expression will be 
used in estimating the thermalization rate of axions in Sec. IIIII 



E. Interaction with other species 



So far we have considered only the self-interaction of the axion field. The highly degenerate axions may also couple 
with other species, such as baryons, relativistic axions, and photons due to the gravitational interactions, though the 
coupling with relativistic axions has already been included in the formalism described in the previous subsection. In 
Ref. [4]. it is claimed that axions have thermal contact with other species after they form a BEG. However, as shown 
in this subsection, there are no such effects at least at the first order in perturbation theory. 

In general, the interaction Hamiltonian with other species b can be written as 

= E iK^i^HbUibj, (67) 

ijkl 
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where A fe lJ fe , is a constant which contains the conservation law of three momenta, and b\ and bi are the operators 
which create and annihilate a species b with the momentum pi, respectively. Here we consider only the interactions 
conserving axion number at the leading order because the interaction rates of all of the axion number violating 
processes are extremely suppressed. It should be also noticed that b particles are conserved as well. 

Substituting Eq. (|67p into Eq. (|54p enables us to calculate how the occupation number of axions evolves with time 
due to the interactions with other particles. The b particles are assumed to be in number states with a distribution 
■A/"&,k(io) at the initial time, 

\{Af},M,W b }) = U ,. r \,„ ■(bl)^ k \Wh{a}), (68) 

k y/Jvb,k'-V b ' k 

where \{Af}, {a}} is given in Eq. ((55)1. A4,fc can take a large number if b is a boson, but it takes either or 1 if b is a 
fermion. In the following, the state (|68[) is taken in computing the expectation value in Eq. (f54j) . Noting that 



jkl 

we find for the condensed modes, 



Hi, b (t),M P ] = E KV-^M&k&i - H.c] , (69) 



([H!, b (t),rt P ]) = iEE [K P \^ t(Ep - Ek)t ^ p al - c.c] for p < K, (70) 

3 k<K 

which exactly vanishes because of the conservation law of three momenta S p+ j,k+j in ^b^kj- This term also vanishes 
for the particle-like modes, 

([H Iti (t) i fi p ])=0 for p>K. (71) 

Hence there are no contributions from interactions with other species. 

From the above discussion, we conclude that the scattering does not occur, in general, between particle-like modes 
in the tree level of the interaction. This is inevitable consequence that follows from the two assumptions, the ^-number 
conservation in Eq. ((67)) and the number state representation for b particles in Eq. ([68)) . Momentum transfer does not 
occur between number states in the tree level because of the conservation of three momenta. Schematically, this fact 
can be understood by using diagrams shown in Fig. [TJ In the usual calculation of S-matrix for the scattering process 
a + b — > a + b, we specify "in" and "out" states as definite particle states for a and b species, but the momentum 
of each particle can differ between "in" and "out" states, as shown in Fig. [T] (a). In the in- in formalism, this tree 
level diagram is deformed such that "in" and "out" states are synchronized [see Fig. (D(b)]. Then, the momenta of 
two "in" states must be the same if they are represented as number states. It is obvious that there is no momentum 
transfer in such a process, and hence it is forbidden. On the other hand, the momenta of "in" states can differ if 
they are coherent states, since they are not eigenstates of the number operator. Therefore, the tree level process is 
allowed for self- interactions between coherent states, as shown in Fig. [1] (c). This is why the second line of Eq. ((62)) 
has a non- vanishing contribution. Figure [T] (d) shows that the scattering between particle-like modes can occur in 
the higher order in perturbation theory. As will be seen in Appendix (A) at least at the second order in perturbation 
theory, this process corresponds to what we calculate by using the usual Boltzmann equation. 



III. FORMATION OF AXION BOSE-EINSTEIN CONDENSATION 



In this section, we discuss the cosmological evolution of dark matter axions based on the results in the previous 
section. The zero mode axions are produced at the time t\ satisfying the condition 

m(T x ) = 3H(h), (72) 

where H(t\) is the Hubble parameter at t\ and T± is the temperature of radiations at that time. The temperature- 
dependence of axion mass is obtained in [2l| , 

A 4 / T \ ~ n 

TO 2 (T) = L68xl0 -7_|CD (73) 

Fi V a qcd/ 
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FIG. 1: Schematics of interaction processes, (a) The usual Feynmann diagram for the tree level scattering process a + b ^ a + b. 
The momentum transfer occurs at the vertex denoted as A. (b) The tree level diagram for the process a + b — s- a + & in the 
in-in formalism. If b is a particle-like mode, two momenta of "in" states at t = to must be same. No momentum transfer 
occurs at the vertex, (c) Tree level diagram for the self-interaction of condensed modes a + a—s-a + a in the in-in formalism. 
Momenta of "in" states can differ from each other. Momentum transfer occurs at the vertex, (d) The diagram for the second 
order scattering process between particle-like modes b + b — > b + b in the in-in formalism. Momenta of "in" states at t = to are 
same, but momentum transfer occurs at each vertex. 



with n = 6.68 and A QCD = 400MeV. The time t\ is estimated as [13] 

t 1= 3.01x10 sec( — ) ^-^) (ir^MeV J ' (74) 

where <7* : i is the radiation degrees of freedom at the time t\. This time scale t\ should be identified with to, which 
was used in the previous section. From Eq. ([50)) , the number density of the zero mode axions is estimated as 

. 2.14 X 1* W* (*) <4+2 '■ ,/<4+ " , (J^Y(«hl)\ m 

V 70 / Vl0 12 GeVy V400MeVy \ R(t) J ' V ' 

where X is a numerical factor determined by the initial misalignment angle (6* lm ) 2 and R(t) is the scale factor of the 
universe at the time t. Since the momentum dispersion of axions is given by the horizon at the QCD phase transition 
5p(ti) ~ , their velocity dispersion is estimated as 

Sp(t) 1 (R(h)\ 4 f g^ Y /2(A + n) / Fa y V( 4 +") ( R{tl )\ 



13 



where we have used the expression for the zero-temperature axion mass given by 

A 4 

m 2 (0) = 1.46 x 1(T 3 ^|^. (77) 
One can easily verify that the state occupation number of the zero mode axions is huge, 

(Itt) 3 rn , \-n/(4+«) / F \2(8+n)/(4+n) / . >. -4 

"-"T^k^^m (lo4v) (iigv) ■ < 78 > 

As was assumed in Sec. Ill C[ if states around the ground state are occupied by AT particles. Therefore, each state 
is occupied by Af/K particles on average. Their thermalization rate is given by the time scale of the change of the 
occupation number for condensed modes Q, 

~ Af p (t) dt ' K 1 

where Af p (t) = (Af p (t)) is given in Eq. ([54| and can be estimated by using the formalism described in the previous 
section. Axions form a BEC if this thermalization rate exceeds the expansion rate H(t) Q. Here, it should be kept in 
mind whether the system is in the condensed regime or in the particle kinetic regime. For this purpose, it is necessary 
to compare T with the typical energy dispersion of axions, 

* ~ i„<o )Wt »* ~ 3, 2 , ur^v m — (*) " _ "" l4+ " , mrm- « 

Let us estimate the thermalization rate in the condensed regime. Substituting the time derivative of Eq. (|64l) into 
Eq. ([79]) yields 

-L condensed — ~T~7~- 77o 

M -Qffca;a*a*] , (81) 

p j,k,l<K 

where we have used the approximation (|66| . Defining the factor A such that 

A^=AVS k+l , p+j , (82) 
and using the approximation J\f p ~ \a p \ 2 — Af/K, we finally obtain 

N 

Tcondcnscd ^ A— = An(f), (83) 

where n(t) is the number density of the zero mode axions given in Eq. ()75j> . For A0 4 type self- interaction, the 
expression of A in Eq. (|56"| gives 

1 condensed, s — ^ 2 ' I / 



On the other hand, the expression of A in Eq. (f58|) for gravitational self-interaction leads to 

_ ^Gm 2 n(t) 

-L condensed, o — / e J \""/ 

m(t)) 2 

where <5p(t) is the momentum dispersion of axions [see Eq. (|76|)]. Note that these expressions are valid only if the 
condition 5ui <C r con( jenscd is satisfied. 

In the opposite case Sco ^> T, the expression in the particle kinetic regime must be used. As seen in Sec. Ill Dt the 
first order term vanishes in this regime, which requires us to evaluate the second order terms in order to estimate 
the transition rate of the occupation number. Since the thermalization rate is the quantity of C(A 2 ), it is suppressed 
compared with that in the condensed regime by a factor of A. Thus, one can expect that it is difficult for axions to 
thcrmalizc in this particle kinetic regime. In fact, it is also obvious from the fact that the conditions r part i c ie > H 
(thermalization condition) and r par ti c io < 5uj (particle kinetic condition) are incompatible due to Slo < H after the 
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time ti [see Eqs. (|T0[) and ([80])]. Here, r par ti c ic is the transition rate obtained from the second order terms in the 
perturbative calculation. Thus, we can conclude that the axion thermalization occurs only in the condensed regime. 

Figure [2] shows the time evolution of thermalization rates T together with the expansion rate H. We find that the 
transition rate r con <3 C nscd,g due to the gravitational self-interaction exceeds the expansion rate when the temperature 
of photons becomes 



Tbec - 2.07 x 10 3 eV X 



V 70 ) 



-3n/4(4+n) 



F n 



10 12 GeV 



6/(4+n) 



( Aqcd \ 
I 400MeV / 



(86) 



which corresponds to the time scale 



t BEC - 3.09 x 10 5 sec X~ 2 (^) 



3n/2(4+n) 



10 12 GeV 



-12/(4+ra) 



( Aqcd V 
I 400MeV / 



(87) 



Before the time £beCj axions are decoupled from each other and are described as a classical field. Once a BEC is 
formed, however, axions behave like cold dark matter by different reason from that of the classical field. In particular, 
from the causality, the correlation length I of the axion field is expected to extend over the horizon I < t . Hence 
the momentum dispersion 8p appearing in Eq. (|85[) becomes comparable with Z _1 ~ t , which makes the time scale 
of the thermalization process much faster. Then, the axion BEC continues to re-thermalize itself, and almost axions 
stay in the lowest energy state. This leads to the modifications of some quantities such as the energy-momentum 
tensor and the evolution equation of density perturbations, but they do not induce any effect on the length scale 
relevant to cosmological observations [!, . 
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FIG. 2: The time evolution of the relaxation rates V in the condensed regime and the expansion rate H . r con d C nscd,s is the 
relaxation rate due to X(f> 4 interaction estimated in Eq. (|84[) and r con d C nscd, 9 is the relaxation rate due to gravitational self- 
interaction estimated in Eq. (|85p . In this figure we normalize all the scales in unit of t\ = 1, where t\ is given in Eq. (|74[1 . and 
the following parameter values X = 1, g*,i = 70, and F a = 10 12 GeV are taken. Axions form a BEC at the time £bec given in 
Eq. (HZ). 

Once the zero mode axions form a BEC, it is expected that they establish the Bose-Einstein distribution with a 
temperature T a different from the photon temperature T. This fact implies that a fraction of axions stay in the 
higher energy states. Then, we need to treat the thermally excited modes and zero modes separately, and the number 
density of axion BEC is given by no -I- ht, where no and ut are the number density of zero modes and excited modes, 
respectively. The total energy density of axion BEC is also written as mno -I- pt with px being the energy density 
of thermally excited modes. Since the number density and energy density of axions just before the formation of the 
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BEC is given by n(iBEc) in Eq. (j75|) and mn(tBEc), respectively, the conservation of number and energy of axions 
yields px = mnx, which shows that the thermally excited modes are non-relativistic. Since the momentum dispersion 
of axion BEC is given by the inverse of the horizon <5j?(£bec) = mSv(tBEc) ~ *bec> the temperature of axions at fsEC 
becomes 



T a (i BEC )~ 3m(Q) -2.48x10 eVI ) [w^eV ) 



V QCD 

400MeV 



which is much smaller than the mass of the axion m ~ 10 6 -10 2 eV. The number density of the thermally excited 
modes is estimated as 

■mT a (t B Ec)\ 3/2 -( m -u)/T„(tT,KC.) ~. ( TOT a(*BEc) 



n T (T a (t BEC )) = ( wT ^ BEc) ) e-(™^)/^(*B E c) 



V 2tt 

, „ v — fln/2C/t-4- 

1.54 x 10 - 50 cm" 3 X 



, a , N -9n/2(4+n) / F \36/(4+m) / . n6 

\70J {lO^GcVj \400MeV) ' 1 ' 

where fi is the chemical potential and we used (m — p)/T a <C 1 for the highly degenerate case. We see that the 
population of thermally excited modes is negligible compared with the total number of axions, 

U T (T a (t BEC )) ^ x lQ _ SOx2 f g^ V 2n/(A + n) f Fa ^(4-n)/(4 + „) ( ^ y 



V 70 / ll0 12 GcV/ UoOMeV/ ' { ' 



n(t B Ec) V 70 J \10 12 GeV J V 4 00MeVy 

Hence, almost all axions are stay in the lowest energy state. 

In Refs. [U, Q, it was pointed out that these excited modes would be relativistic and affect the cosmological 
parameters such as the baryon to photon ratio and the effective number of neutrino species, if axions enter into 
thermal contact with photons. However, as shown in Sec. Ill El the interaction between axions and other species 
exactly vanishes at the first order in perturbation theory, and hence the thermalization rate with other species is 
heavily suppressed, which indicates that axions forming a BEC are decoupled from other particles and do not give 
any significant modifications to cosmological parameters. Even though, if axion BEC is the dominant component of 
dark matter, they give some imprints on the structure of the inner caustics of galactic halos 17']. This might be a 
useful tool to distinguish axion dark matter from other particle dark matter candidates. 



IV. SUMMARY AND CONCLUSION 



In this paper, we developed the formalism to describe the evolution of the zero mode axions in terms of the quantum 
field theoretic method. In order to solve the evolution of occupation number of axions, we used the in-in formalism 
and the coherent state representation for a highly degenerate Bose gas of axions. Combining these two ingredients, 
we derive the time evolution of the expectation value of the number operator in a perturbative way. We showed that 
there is non- vanishing contribution for self-interaction of condensed modes at the leading order in perturbation theory 
[see Eq. (|64p j. On the other hand, the interactions between particle- like modes including other species exactly vanish 
at the first order in perturbative expansion, which indicates that their relaxation rate is suppressed. 

Using the results for the time evolution of the occupation number, we estimated the thermalization rate of the 
zero mode axions. We recovered the expressions for thermalization rates obtained in and confirmed that axions 
form a BEC due to the gravitational self-interactions when the temperature of photons becomes O(10 3 )eV. After 
the formation of axion BEC, almost all axions are in the lowest energy state and they continue to re-thermalize 
themselves. 

From the fact that the tree level contribution vanishes for the interaction with other particle-like species, we conclude 
that the axion BEC has no thermal contact with other cosmological fluids, and that there is no significant effect on 
cosmological parameters. In particular, for the effective number of neutrino species N e s, the result of [4j, |8| predicts 
the higher value N c g = 6.77 than the standard model, but in our analysis N c s does not differ from the standard value. 
Hence the axion BEC is consistent with the standard cosmology. Only peculiar prediction is the specific phase space 
structure for galactic halos |7|, which gives a possibility to probe axion BEC dark matter on observational grounds. 

Finally, let us comment on a speculative point in our discussion. For the gravitational self-interaction of axions, 
we use the expression (|57|) which holds in the Newtonian limit. The use of this term is justified only if we are able 
to ignore the requirement of the causality. To be more precise, Eq. (1571) can be regarded as a good approximation 
while the time scale of the interaction T _1 exceeds the typical length scale SI ~ (5p _1 on which the interaction takes 
place. After the formation of axion BEC, however, it is expected that Sp ~ H, and hence this condition becomes 
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r _1 > H" 1 , which seems to be incompatible with the thermalization condition V > H. Therefore, in order to 
make a clear description about the rethermalization of axion BEC, we must extend our formalism into the expanding 
background, including the correction coming from general relativity. This issue is left as our future work. 
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Appendix A: Second order in perturbation theory 



In this appendix, we compute the terms of second order in Hj [the third term of the right hand side of Eq. (I54[) ] 
From Eqs. (|55l) and (|60|). we obtain 



[Hi(t),J\F p ] = -L [^e^alala.a, - H.c. 



(Al) 



jkl 



Then we find 



[ff/(ti),[ff/(t2),Ay 



i i 

8 V* 



£ A-»A^ e -^r^+^^) [a t a t amania t a t a . ap] 



mnqrjkl 



l l 



SVB E A™A^.e-^- ^- n l^[alaia m a n ,alata kai } 



(A2) 



rnnqrjkl 

Note that 

[a\ala m an, a\a\aja p ] = V 2 (S n iS mk + 5 n k8 m i)a\alaja p - (S jq S pr + 5 P q5j r )a\a\a m a n 



-V 



j j t | x T T T X ^ T T X^^t 



t t t 

7 I rt 1 r\ ' 



t t t 

1 rt ' n i t 



t t t 

i ' rt ' rt < i 



+8 n ia { q ala\a m a J ap + 8 n ka q alaJa m a 3 a P — S jr a q al.a]a m apa n — 8 pr a q a\a}a m aja n . (A3) 
Using this formula, after some simplification, the expectation value of Eq. (|A2|) reduces to 



t t t 



t t t 
I rt 1 n f 



t t t 

i • rt n ' f 



([HjihUHzih),^]}) 



1 1 
1 1 



qjklm 



qjklm 



1 l 

2 V? 
1 1 



£ [A^e-^^"^ (alalala p a 3 a n ) 



nmqjkl 



i^f E [ArAS e - i( ^" il+ < t2) (4atata p a m a n )+c.c 



7imqjkl 



l l 

~4 



]T [A™A^ e -(^T* 1 + f2 "^)( a t a t a t aj . aman) + c . c .] . (A4) 



nmqjkl 
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To compute the expectation value, we note that 



a k >>a k ,a k \{Af}, {a}} 



^/AfkWk - 1)(A4 - 2)^ 3 / 2 |{A0 3fc , {a}) 



y/N k N k >{N k > - l)V*/ 2 \{N} k ' 2k ' ,{a}) 
^N k N k ,N k „ V^ 2 1 {AT}"*' < fc " ,{a}) 



- l)a k „VV 2 \{N} 2k ,{a}} 
V7JJ^a k „V 3 / 2 \{M} kk ',{a}} 
^W k a kl a k „V^ 2 \{M} k ,{a}) 
a k a k ,a k „VV 2 \{N},{a}} 



f k = k' = k" > K 

f k ^ k', k' = k" > K, and k > K 

f k^k' ^ k", k> K, k' > K, and k" > K 

f fc = k' > K and k" < K 

f fe ^ k', k > K, k' > K, and fc" < if 

f k > K, k' < K, and fc" < K 

i k<K, k 1 <K, and fc" < K 



where the state |{A/'} 3 ' C , {a}) contains 



(A5) 

factor (4) JV '' i " 3 /v / P4 - 3)!^- 3 for mode fc, the 



state |{A/"} fc < 2fc ',{a}) contains a factor (4) M, ~ 1 (°fc')' A '*'~ 2 /'^C^b - l)!^*- 1 ^' - 2)!V^*'- 2 
for modes fc and fc', and the state |{A/"} fc ' fc ,fc , {a}) contains a factor 
(al)**- 1 (al)** , - 1 (al l ) Jsf k"- 1 /y/(Afk - l)!^*- 1 ^-' - l)!^*'- 1 ^ - i)|yJV fc »-i for modes fc, fc' and fc". 
Using Eq. (|A5|) . the expectation value of Eq. ()A4[) can be evaluated in a similar way to the first order terms. With a 
tedious but straightforward calculation, the first line of Eq. (|A4[) becomes 



1 1 

lye 



c.c. 



E [ArA»e-^ fl+ ^ t2 ) (a t fl t ajflp>+ , 

qjklm 

' ' - \ \ \ I Y!\" \ 

qm 

qm k<K 1>K " 

1 1 E E [ a ^a- e-^ir^Lv alatajap 



4 V 4 



?ra j, k ,l<K 



for p < K, 



(A6) 



1 1 
4V^ 



qjklm 

l l 

It74 



E [ArA^e-^r* 1+ n- t3)(a t a t a . ap>+c . 
E [A^A» e-W^+nSS.*"^^ - 1) + c. 



-5^1 E E l A p?l 



9" 3>K,j=£p 



1 1 

2^ 



gm j, k <K 

The second line of Eq. (|A4[) becomes 



1 1 



E ["p'" 

qjklm 
1 1 



for p > K. 



A^ ro A^e-<^+^r^)(at a t a . a9)+c . 

AWA^e-'^+'C t2) AA 9 (AA 9 - 1) + . 



e- m ^-*'^^ q + c.c. 



iivy [ 

rn q>K 

2^4 E/ E/ |Apml 

m j>K,j^p q>K 

•4E E E [A^AS-e-^'^^'^a^M 

m j,k<Kl>K 

i^E E [A^ASTe-^^^afaja. + cc. 

m j,k,l,q<K 



(A7) 



(A8) 



The third line of Eq. (|A4[) becomes 

\jf E [A^A^e- < ^* 1 +^^)( a t a t ia t apa . an> + c . c . 



or 



nmqjkl 

1 

v 1 



!iEEE [A;%e-««"<^^(^-i) Q ; ap+ c. 

Z g <if fc>A' 

^EEE [Ai fe fe A^e-^V I+ n- t3 )^ ( ^_ 1)a * ap + c 



Z q<Kk>K 

71 EE E E [A^A4a*a p 



F 4 



( q<K m=£k,m>K k>K 



( A ik A pm -i(Q' g k m tx+Q p J-t 2 ) , A ;m A pfe p-icn^ti+njl-ta) , A ifc a p™ p-«(^ fe ti+np"*t 2 )\ , 

1 gm n fc| L ' ^ 1 ^qm ± ^kl ' 1 ^mk li -ql I 1 



+^E E E [^<^a p a„(A^ £ 



/ m,7i,k<K q>K 



+1^ E E E K<^a p a„ (Ajfc.A^^^) + A^A^ 1 * 



c.c. 



/ m,n,k<K q>K 



e" < Km*i+ n k!*2) a * a ^ a * Q , j)aj .Q !n + c.c. 



for p < K, 



l q,m,k,j,n<K 



E ^Afje- i (* fl ^ il '(atat fl t ap¥n ) + c.c, 

nmqjkl 

= \h E [AfeAjye-^-^^C^ - 1)(AT, - 2) + . 



qjtp,q>K 



1 1 



q^p,q>K 



•^E E [^e-^^^,(V 4 -l)Ar p + . 

/ q^p,q>K 

1 1 

"2 v~ 4 



^E E [A^A^e-^!S^*-W,(^-l)Ar p + c.c. 



1 



EE E [lA^pe-^^-^AA^AA.+e.e. 



/ q^p,q>K k^q,k^p,k>K 



+ 171 E E E A^e^^^WpA^+ac 



; q^p,q>K k^q,k^p,k>K 



V 4 

t/iE E E [a^^^^^)w+c-c. 

( q^p,q>K k^q,k^p,k>K 



^4 E E " lK«m (A^A-re- 4 ^^) + A£A£ e -«**W) + c. 

! n.m<K 

E [AW-l)a> m (A;rX7^^^^^ 



1 1 

"2F 4 



/ n,m<K 



^E E E ^ 



g a n o: m 



n,m<K q^p,q> K 
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1 "n? pi'' 1 p? til 1 pq nl J 



i n.m,j,k<K 

l l 



The fourth line of Eq. (|A4|) becomes 
1 1 



for p > K. 



i^r E [ArAge^ n Sr*^^<4atata p a m a B > + c.c. 



nmqjkl 

i i 
i i 



£ [A^A^, e -^ fc ^^^^(AA fe -l)a; ap + c.c. 



i q<Kk>K 



ITTlEE E [A^Afie-<^^+^^)A4(AA fe -l)a*a p + 

i q<K k>K 

^4 EE E E [AS?Aje-^« ^«^ fe A/>^ + , 



i q<K m^k,m>K k>K 



^EE E E [A^A^e^^*-)^ 



'k-N m a q oip + c.c. 



I q<K m^k,m>K k>K 



^E E E [ A ^A^„e~ l(n "'* 1+n « ; " t2) A^a I *„<a p afe + . 



5tt?E E E[ A ^ e " ,(<tl+Ct2) ^>:^+' 



H m,n,k<K q>K 



I m,n,k<K q>K 



i^I E E [ArAge- i(n 2"^ +n S*^a^*a*a p a m a„ + c.c. 



/ q,m,k,j,n<K 



(A10) 



for p < if, (All) 



1 1 

'4 



£ [A3»Aj5e-*< n 5-*^*«)<aio5aJfl p a m a B > + c.c 
nmqjkl 

E [AfeAjJe-^na^^CArp - W - 2) + c.c. 



5^4 E E [A^e-^?^^^(A^-l)A^ + , 
i q^p,q>K 

^E E [A%e-«'H(Afp-lK 

! q^p,q>K 

1 1 



171 E E [A p ?A^ e -^?^<^)^(^-l)A^ + 1 

i q^p,q>K 

1 1 

2T 3 



-^E E [A^e-^^+<^AA g (AA g -l)AA p + 1 



Z q^p,q>K 



HE E E [AKe-^^^i^M 



V 4 
1 1 



fc + c.c. 



I q=£p,q>K k^q,k=£p,k>K 



E E [|A^| 2 e-<^-^A^A^ + c.c. 



2 V 4 
1 

'V4 



i q^p,q>K k^q,kytp,k>K 

^E E [A;Xe- l( ^ fl < t2) ^(^-l)a; Qm + l 

Z n,m<K 
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i n,m<K 

^lE E E E^pA^c 

/ n,m<K q^p,q> K 



l l 

"2l74 

l l 

Iv 7 



/ n,m,j,k<K 

Finally, the fifth line of Eq. ()A4j) becomes 
1 1 

"If 7 



E E [Ai?A^e" i(n - ^^A^a^ + c. " 

Z n,m,j : k<K 

E E [^^^^-^^^0,0,,+ c. 



for p > K. 



(A12) 



I^r E [A^Age^^+ n fi*"><ajtotata J a m a n >+c.c. 



nmqjkl 

i i 



4 V4 E [A"A|Je- < ( n K^+ n ;S t 'W,(^ - 1)(AA 9 - 2) + c.c. 



q>K 



^lE E [AgA-e-^^^^AA g (AA g -l)AA fe + < 



q>K k^q,k>K 



1 1 



, r * E E [A^A-e-^<^^^)^(^ - 1)A4 + c.c. 



g>if k^q,k>K 



<?>if k^q,k>K 



-A E E [Ar fc A^e- l(n - tl+ °" t2) AA,(AA 9 - 1)A4 + c.c. 



<j>if k^q,k>K 



1 



EE E fA^A-e-^^^^^V^ m + c.c. 



g>i-C k^q,k>K m^k,m^q,m>K 



" 2 y 4 E E 



5] [|A£™| 2 e~ in - {tl - ta) Af q Af k Af m + c.c. 



q>K k^q,k>K m^k,m^q,m>K 



!jE E [A™^e-W^+^^(^-l)<a m + c. 



q>K n,m<K 

l l 

'2F 4 



5^4 E E [AW - (Ag^e-*^^*-) + A^e-^^W) + c.c/ 



q>K n,m<K 



l l 
'Iv 4 " 



E E [A« A^e-^^-^*-)^, - lRa m + c.c. 

q>K n,m<K 

E E KA4a> m (Af fc A^ 



^ e -»(n|*ti+ns™t2) 



g>_fC n,m<K k=£q,k>K 

. \mkApq p -t(n™ fc fc ii+nP«i 2 ) , A ,nk A Pq i^tx+f}^) , Am* AM „-*(n™ fc ti+n^t 2 )\ , 
1 • Ii pfc ls -nq^ 1 pq nk 1 pn ^''qk J 1 

5^4 E E [Af„A-e- l ^ tl+ ^ t2 V^ a > m + c.c. 

q>K n : m<.K 

^E E [Ay^W^*-)^afa m a B + . 

q>K k,l,m,n<K 

\h E E [^«i«r«m«» (AJTAge-W^^) + AJSAre-*^"^') + c.c. 

g>iC k,l,m,n<K 
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\h E E [A^^^-^^aiof a^a,. + , 



1 1 

i 

q>K k,l,m,n<K 



-\yi E [A™Xfe- i(n ™" tl+nS " t2) «^ra*^a m a n + c.c.] . (A13) 

In the particle kinetic regime, these tremendously long equations can be simplified as follows. For p > K, we obtain 

M p {t) ~ yV p (to)- / dt 2 [' dt^lHifallHifa),^]]). (A14) 

J to ./to 

In general, terms which contribute to the expectation value take a form 

Taking the time derivative after performing the integration over t\ and we find 



dJ\f p 



_L e -i(f!i+n 2 )i ^_ e -i(nit +o 2 t) _|_ c c (A15) 



If 7^ 0, these terms drop out because of the rapidly oscillating factor in the particle kinetic regime (£lp l q t —> oo). 

On the other hand, if f2i + f22 = 0, the first term of the right hand side of Eq. (|A15[) cancels with its complex conjugate. 
Then we obtain 

dM p 2 . n . 
—j— oc — sinS2i(t - to), 
at 1 1\ 

Note that the energy conservation emerges in the limit because ft%(t — t ) — > oo, 

-^-sinOi(t-io) 2n8{Cl 1 ), (A16) 
Hi 

which implies that terms with Qi ^ do not contribute to the final result in this limit. For example, the second line 
of Eq. (|A7j) gives Qi = fi£™, which does not vanish because of the conservation law of three momenta in A^™. The 
exception is the case with q = m = p, but the careful inspection shows that this term exactly cancels with the second 
line of Eq. (|A8[) . Similar discussions are applied for the second, third, fourth, and fifth lines of Eq. (IA10[) . the second 
and fourth lines of Eq. (|A12j) . and the second, third, and sixth lines of Eq. (|A13|) . After all, the remaining terms lead 
to 



-L ]T |A^| 2 27n5(^) [A/^ (A17) 



dMp 

~~dT _ 2V A 

klq>K 



where we have neglected the contribution that contains the integration over condensed modes (i.e. ^2 q<K X)fc i>k)' 
because such terms are prohibited by the energy conservation (|A16[) . In this way, we recover the usual Boltzmann 
equation jj]. 
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